Past, present, and future predictions on the suitable habitat of the Slender racer (Orientocoluber spinalis) using species distribution models

Abstract Species distribution models (SDMs) across past, present, and future timelines provide insights into the current distribution of these species and their reaction to climate change. Specifically, if a species is threatened or not well‐known, the information may be critical to understand that species. In this study, we computed SDMs for Orientocoluber spinalis, a monotypic snake genus found in central and northeast Asia, across the past (last interglacial, last glacial maximum, and mid‐Holocene), present, and future (2070s). The goal of the study was to understand the shifts in distribution across time, and the climatic factors primarily affecting the distribution of the species. We found the suitable habitat of O. spinalis to be persistently located in cold‐dry winter and hot summer climatic areas where annual mean temperature, isothermality, and annual mean precipitation were important for suitable habitat conditions. Since the last glacial maximum, the suitable habitat of the species has consistently shifted northward. Despite the increase in suitable habitat, the rapid alterations in weather regimes because of climate change in the near future are likely to greatly threaten the southern populations of O. spinalis, especially in South Korea and China. To cope with such potential future threats, understanding the ecological requirements of the species and developing conservation plans are urgently needed.

Over the last 8000 years, the average annual temperature in Asia rose by about 1°C, but 0.6°C of this increase happened during the last 100 years, showing an acceleration in the pattern (Ennis & Marcus, 1996;Hughes, 2000;Walther et al., 2002). As a result of climate change, species abundance change (Ehrlén & Morris, 2015;Parmesan & Yohe, 2003;Root et al., 2003;Thomas et al., 2004), distribution range shift (Bellard et al., 2012;Habibzadeh et al., 2021), and extinction rate of populations (Bestion et al., 2015;Thomas et al., 2004) are greatly impacted and might strongly fluctuate. Climate change is expected to further accelerate (Perkins et al., 2012), resulting in even more severe impacts on numerous species (Pearce et al., 1996).
Reptiles are one of the animal groups sensitive to environmental alterations because of their ectothermic characteristics and low dispersal ability (González-Fernández et al., 2018;Huey, 1982), and they consequently face greater threats (Sinervo et al., 2010;Urban, 2015). Considering the worldwide reptile decline following the acceleration in climate change (Reading et al., 2010;Winne et al., 2007), developing conservation plans to prevent the decline of species is necessary. A step in this direction, and in particular for reptile groups with small population size, is to conduct basic ecological research such as determining species' distribution and the underlying environmental factors. In addition, predicting the species response to climate change is also essential to understand potential threats to the species and developing long-term conservation plans.
In addition, research on the adaptation of species to past climate change provides a significant perspective on how species will react to climate change in the future (Pearson, 2006). thus becoming an increasingly important tool for ecology, evolution biology, conservation biology, and climate change biology (Botkin et al., 2007;Guisan et al., 2013;Mi et al., 2017;Zhang et al., 2011).
Specifically, these modeling approaches across past, present, and future can be most important to understand the status of threatened, rare, and less known species (Araújo et al., 2004;Fois et al., 2018;Sinclair et al., 2010).
The Slender racer (Orientocoluber spinalis) is the only species within the genus Orientocoluber. The species was recently transferred from the genus Hierophis to Orientocoluber based on morphological, osteological, and biogeographical characteristics (Kharin, 2011;Nagy et al., 2004). This species has a wide distribution through cen-  Kharin & Akulenko, 2008;Munkhbayar & Munkhbaatar, 2012). Field observations of the species are very rare across all countries (Kharin & Akulenko, 2008;Kim & Han, 2009) (Anderson & Raza, 2010;Radosavljevic & Anderson, 2014;Kim, Park, Bae, et al., 2020) using a 1 km, 60 km, 75 km, and 100 km radius. For each radius, we selected only one random presence datapoint among the overlapping data. We calculated the density of filtered presence datapoints and visualized their distribution with the "heatmap" tool in QGIS (QGIS.org, 2020). The distribution of the presence data was the most appropriate based on the 75 km radius when compared with the distribution of the species. Thus, for further analyses, we retained a total of 94 presence datapoints (16 in South Korea, 1 in North Korea, 72 in China, 3 in Mongolia, 1 in Russia, and 1 in Kazakhstan; Figure 1, Table 1).  We performed the SDMs for past and future prediction based on the Community Climate System Model 4.0 (CCSM4), which is often used for amphibians and reptile SDMs (Borzée et al., 2019;González-Fernández et al., 2018;Kim, Park, Bae, et al., 2020). For the past modeling, we used three time windows: (1) the last interglacial period about 130,000 years ago, when the average air temperature was 5°C higher and sea level was 2.2 to 3.4 m higher than present (Andersen et al., 2004;Otto-Bliesner et al., 2006). were similar to the present (Gagan et al., 1998). SDMs based on these historical climatic predictions have been conducted in various animal groups, including G. japonicus ,

| Modeling methods
We performed the SDMs on O. spinalis using a maximum entropy modeling software (MaxEnt v. 3.4.4;Philips et al., 2020) because of its high predictive power in both extrapolation and interpolation with presence-only data (Heikkinen et al., 2012;Phillips et al., 2006;Wisz et al., 2008). MaxEnt is one of the most efficient approaches for predicting predicting species' potential distributions (Elith et al., 2006;Elith et al., 2011;Phillips et al., 2006). We ran the models with 25% random test, 15 bootstrap replicates, 5000 iterations, one regularization multiplier, and a logistic output (Phillips, 2005;Young et al., 2011). We used the area under the curve (AUC) and the true skill statistic (TSS) to validate the model reliability using the sdm package in R (Naimi & Araujo, 2016). The AUC is based on the receiver operating curve (ROC), and it is commonly used to evaluate model performance (Allouche et al., 2006;Bradley, 1997). TSS is a simple and easy method to verify the sensitivity and specificity predictive power of SDMs (Allouche et al., 2006;Shabani et al., 2018).
To determine the suitable habitat and its rate of change for O. spinalis, we designated a threshold using the approach called "maximizing the sum of sensitivity and specificity (max SSS)" as it is an adequate method for presence-only data (Liu et al., 2013). In addition, we used the jackknife test and computed the response curves built into the MaxEnt software to evaluate the contributions of each environmental variable in the models (Jiménez-Valverde, 2012). Also, to identify climate extrapolation in different periods, we used the multivariate environmental similarity surfaces (MESS) analysis incorporated into MaxEnt software (Archis et al., 2018;Elith et al., 2010). As our past and future models were built based on present climate variables , it is necessary to be careful with interpretations for the area outside of the present climate range (Carneiro et al., 2016). and annual precipitation of around 600 mm were critical for the suitable habitat of O. spinalis (Figure 2). The contribution of each variable to the model was in the order of annual mean temperature (49.8%), isothermality (16%), annual precipitation (14%), precipitation seasonality (10.4%), altitude (6.6%), and mean diurnal range (3.3%).

| Climate similarity
The modeled climate for 6 kyr ago and 2070 did not generally deviate from the present climate range; however, several regions were dissimilar 130 and 22 kyr ago, and warrant caution when interpreting the data (Figure 3).

| Distribution shift
According to the present model results

| DISCUSS ION
Here, we show that O. spinalis has a wide distribution, unlike related species that occur across a relatively restricted range (Das et al., 2019;Mirza et al., 2016;Nagy et al., 2004). For example, Hierophis spp. mainly inhabit southern Europe, a warm Mediterranean climate (Jablonski et al., 2017;Rato et al., 2009) (Kottek et al., 2006;Peel et al., 2007), with an annual mean temperature of around 9°C and annual precipitation of around 700 mm. The continental climate with high habitat suitability for the species is also characterized by a large annual temperature variation (Kottek et al., 2006;Peel et al., 2007). Therefore, appropriate isothermality could be another  (Peel et al., 2007). In opposition, the high-altitude mountain areas and Considering that O. spinalis is not currently found in Japan and con- in the distribution of many species (Hughes, 2003;McLeman & Smit, 2006;Pearson, 2006), and more than 68% of species are likely to see their habitat shift poleward (Karl et al., 1996;Thomas, 2010;Walther et al., 2002). Regionally, mid-southern areas of the Korean Peninsula will lose suitable habitat for O. spinalis by 2070. SDMs studies on species distributed in South Korea such as G. japonicus (Kim, Park, Bae, et al., 2020), Karsenia koreana (Borzée et al., 2019), and Onychodactylus koreanus (Shin et al., 2021) also showed that the suitable habitat will be shrunk gradually because of climate change.

ACK N OWLED G M ENTS
We thank Il-Hoon Kim, Woo-Jin Choi, Dae-in Kim, and Irina V.
Maslova for assistance in data collection. This work was supported by the Basic Science Research Program through the National Research Foundation of Korea funded by the Ministry of Education (2020R1I1A3051885).

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The input data and output data in this study are accessible at Dryad